require(data.table)
library("dplyr")
library(ggplot2)
library(forecast)
##############################
############################  Figure 3a-d time series effects
##############################

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"



#================load data=================
#load dif.8crop.wm prepared by Data.Fig3-4.R
load(paste0(wd1,"/data/RData/Fig3-4.globalmean-yieldchanges.Rdata")) #including all grids with crop failures to calc modeled total change in yield under different scenarios (area>1000ha)
dif.8crop.wm[Year==2090,sum(area), by=.(Scenario,Year)]
dif.8crop.wm[Year==2090,sum(area2), by=.(Scenario,Year)]

#=smooth data (do wtmean on log effects and then smooth over time; similar results as wtmean on smoothed grid data dif.8crop.ma)
dif.8crop.wm.ma <- dif.8crop.wm[, lapply(.SD, ma, order=5), by =.(Scenario,Crop,Irrig)] 
dif.8crop.wm.ma <- dif.8crop.wm.ma[, lapply(.SD, as.numeric), by =.(Scenario,Crop,Irrig)] 
dif.8crop.wm.ma[,Year:=as.integer(Year)]
dif.8crop.wm.ma <- dif.8crop.wm.ma[!is.na(Year)]


#load MLR results with confidence intervals from Bootstrap resampling (by Data3.MLR-log-bootstrap-resampling.R)
load(paste0(wd1,"/data/RData/MLR-pergrid-log-wtm-effect-Bootstrap-1000.Rdata"))

#========Start plot time trends of estimated partial and total effects
#=first translate log to % effects and then plot the original or smoothed time series
eff.boots.pct <- eff.log.boots[dif.8crop.wm,
                                    .( Level, T.eff=(exp(T.log)-1)*100, #R.eff=(exp(R.log)-1)*100, M.eff=(exp(M.log)-1)*100,
                                       RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log)-1)*100, P.eff=(exp(P.log)-1)*100, RH.eff=(exp(RH.log)-1)*100,
                                       Tot.eff=(exp(Tot.log)-1)*100, co2.eff=(exp(i.co2.log)-1)*100, luc.eff=(exp(i.luc.log)-1)*100,  
                                       Y.mod.dif=(exp(i.Y.mod.dif)-1)*100, #Y.mod.dif=(yield-yield0)/yield0*100, #differene of wtm yield is not the same as wtm of per-grid difference
                                       Tot.eff2=(exp(Tot.log2)-1)*100, Y.mod.dif2=(exp(i.Y.mod.dif2)-1)*100, #climate+CO2 for ER
                                       Tot.eff3=(exp(Tot.log3)-1)*100, Y.mod.dif3=(exp(i.Y.mod.dif3)-1)*100, #climate+co2+LUC
                                       #Tot.eff3=(exp(Tot.log+i.co2.log+i.luc.log)-1)*100, Y.mod.dif3=(exp(Y.mod.dif+i.co2.log+i.luc.log)-1)*100,
                                       yield=i.yield, yield1=i.yield1, yield0=i.yield0, area=i.area , yield2=i.yield2, area2=i.area2 #
                                    ), by=.EACHI, #on=.(Scenario,Year)]
                                    on=.(Scenario,Crop,Irrig,Year)]
eff.boots.pct <- eff.boots.pct[!is.na(Level)]
#aggregate per crop effect to global considering base yield difference between crop types
eff.tot <- rbind(eff.boots.pct[Scenario=="RCP45 - RCP85",
                                     .(T.eff=weighted.mean(T.eff, yield1*area, na.rm = T), #R.eff=weighted.mean(R.eff, yield1*area, na.rm = T), M.eff=weighted.mean(M.eff, yield1*area, na.rm = T),
                                       RD.eff=weighted.mean(RD.eff, yield1*area, na.rm = T),RI.eff=weighted.mean(RI.eff, yield1*area, na.rm = T),
                                       P.eff=weighted.mean(P.eff, yield1*area, na.rm = T), RH.eff=weighted.mean(RH.eff, yield1*area, na.rm = T),
                                       Tot.eff=weighted.mean(Tot.eff, yield1*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield1*area, na.rm = T), 
                                       co2.eff=weighted.mean(co2.eff, yield0*area, na.rm = T), luc.eff=weighted.mean(luc.eff, yield*area, na.rm = T), #LUC has differnet base yield from RCP45_85LU
                                       Tot.eff2=weighted.mean(Tot.eff2, yield0*area, na.rm = T), Y.mod.dif2=weighted.mean(Y.mod.dif2, yield0*area, na.rm = T),
                                       Tot.eff3=weighted.mean(Tot.eff3, yield0*area, na.rm = T), Y.mod.dif3=weighted.mean(Y.mod.dif3, yield0*area, na.rm = T), 
                                       yield1=weighted.mean(yield1, area), yield=weighted.mean(yield, area), yield2=weighted.mean(yield2, area2), area2=sum(area2, na.rm=T),
                                       yield0=weighted.mean(yield0, area), area=sum(area, na.rm=T)),
                                     by =.(Scenario,Level,Year)],
                  eff.boots.pct[Scenario!="RCP45 - RCP85",
                                     .(T.eff=weighted.mean(T.eff, yield0*area, na.rm = T), #R.eff=weighted.mean(R.eff, yield0*area, na.rm = T),M.eff=weighted.mean(M.eff, yield0*area, na.rm = T),
                                       RD.eff=weighted.mean(RD.eff, yield0*area, na.rm = T),RI.eff=weighted.mean(RI.eff, yield0*area, na.rm = T),
                                       P.eff=weighted.mean(P.eff, yield0*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield0*area, na.rm = T),
                                       Tot.eff=weighted.mean(Tot.eff, yield0*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield0*area, na.rm = T),
                                       yield0=weighted.mean(yield0, area), yield=weighted.mean(yield, area),
                                       area=sum(area, na.rm=T)),
                                     by =.(Scenario,Level,Year)], fill=TRUE)

#==5 year runnning mean of global effects
eff.tot.ma <- eff.tot[, lapply(.SD, ma, order=5), by =.(Scenario,Level)] #by =.(Scenario,Crop,Irrig)] 
eff.tot.ma <- eff.tot.ma[, lapply(.SD, as.numeric), by =.(Scenario,Level)] #by =.(Scenario,Crop,Irrig)] 
eff.tot.ma[,Year:=as.integer(Year)]
eff.tot.ma <- eff.tot.ma[!is.na(Year)]



#=======Linegraph for trends of effects in 21st century

#=======Final Style: partial effects of T,RD,RI,P,RH,CO2 (assume ER is independent of LUC)
#=======aggregate log effects (weighted by area only) before translating to % effect give best results on total and partial effects (see eff.mc.wm.crop)
#=======Do not use pct wm effects or resampling on pct effects, they will bias partial effects of RD and RI although total R.eff look ok, the sum of T.eff+R.eff+M.eff is also larger than Tot.eff!
eff.tot1 <- copy(eff.tot.ma)
eff.tot1 <- eff.tot1[Scenario=="RCP45 - RCP85", Tot.eff:=Tot.eff2] #climate+CO2 only
eff.tot1 <- eff.tot1[, Y.mod.dif:=(yield-yield0)/yield0*100] #relative change of global average yield (yield=RCP45_85LU,SAI,MSB,CCT, yield0=RCP85)

#compared tot.eff and Y.mod.dif
eff.tot1[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff), by=.(Scenario)]
eff.tot1[Year >=2075 & Year<=2099 & Level=="Mean", mean(co2.eff), by=.(Scenario)]
eff.tot1[Year >=2075 & Year<=2099 & Level=="Mean", mean(Y.mod.dif), by=.(Scenario)] 
#statistical estimated Tot.eff (weighted mean of grid level change) is in-between of two averaging methods for CLM5 modeled yield change
#this is 5-year running mean smoothed

#% effect on global yield
eff.tot1.t <- eff.tot1[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Year,Scenario)]
names(eff.tot1.t)[3:6] <- c("Variable", "Mean", "Conf.L", "Conf.U")



#======================select distinct color-blind colors for scenarios, variables, and crops
library(RColorBrewer)
cbPalette <- brewer.pal(n =12,"Paired")
grays <- c("#D9D9D9", "#BDBDBD", "#969696", "#737373", "#525252", "#252525", "#000000")
barplot(height=rep(1,length(cbPalette)),col=cbPalette) #display colors
varcolor <- c(cbPalette[c(10,7,8,1,2)],"#E64B35FF", "#000000", "#737373")
barplot(height=rep(1,length(varcolor)),col=varcolor) #display colors

#scientific journal colors
library("ggsci")
cbPalette <- pal_npg("nrc", alpha = 1)(10) #nature journal colors
# cbPalette <- pal_aaas()(10) #AAAS journal colors
# cbPalette <- pal_lancet()(9) #lancet colors
barplot(height=rep(1,length(cbPalette)),col=cbPalette) #display colors


#====final colors for variables and scenarios
varcolor <- c("#E64B35FF","#FF7F00" ,"#FDBF6F", "#A6CEE3", "#1F78B4", "#00A087FF", "#000000", "#737373")
barplot(height=rep(1,length(varcolor)),col=varcolor) #display colors

#casecolor <- c("#7E6148FF", "#B09C85FF", "#CC79A7", "#8491B4FF", "#91D1C2FF")
casecolor <- c("#440154FF", "#7E6148FF", "#CC79A7", "#8491B4FF", "#91D1C2FF")
barplot(height=rep(1,length(casecolor)),col=casecolor) #display colors


linetype <- c(1,1,1, 1,1,1,1, 1)
size <- c(1,1,1, 1, 1, 1, 1, 1)*1
override.linetype <- 1
override.size <- 1.1
casename = c("'SAI - RCP85'", "'MSB - RCP85'", "'CCT - RCP85'", "'Emissions reduction'")
# casename = c(expression('('*bold(a)*') SAI - RCP85'), expression('('*bold(b)*') MSB - RCP85'), expression('('*bold(c)*') CCT - RCP85'), 
#              expression('('*bold(d)*') Emissions reduction')) #expression('('*bold(i)*') RCP45 - RCP85'))
scenario = c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85")
#verifying the match between regression estimate and modeled change
vars = c("T.eff", "RD.eff","RI.eff", "P.eff", "RH.eff", "co2.eff", "Tot.eff", "Y.mod.dif")
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             #expression(atop(atop(textstyle('R'), textstyle('effect')),' ')),
             #expression(atop(atop(textstyle('M'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('P'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('RH'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             #expression(atop(atop(textstyle('LUC'), textstyle('effect')), ' ')),
             #expression(atop(atop(textstyle('Total'), textstyle('effect')), ' '))) # textstyle tells not to shrink the text size when atop more than 2 layers
             expression(atop(atop(textstyle('Total'), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(CLM5)')), ' ')))

p0 <-  eff.tot1.t[Variable%in%vars] %>%
  mutate(Var = factor(Variable, levels=vars)) %>%
  mutate(Scenario = factor(Scenario, levels=scenario, labels=casename)) %>%
  ggplot(aes(Year, Mean, group=Variable)) +
  #ggplot(aes(Year, Median, group=Var)) + #Weighted median shows better match with modeled yield change than weighted mean
  geom_line(aes(color=Var, linetype=Var, size=Var)) + #use Case as group and set override values in manual scale
  geom_hline(aes(yintercept=0), color="gray") +
  geom_linerange(aes(x=Year, ymin=Conf.L, ymax=Conf.U, color=Var), alpha=0.3, show.legend=FALSE) +
  #change linetype, color, size mannually
  scale_linetype_manual(values=linetype, guide = FALSE)+
  scale_size_manual(values=size, guide = FALSE) +
  scale_colour_manual(name="", labels=varname, values=varcolor)+ # cbPalette[c(1:4,4,5,5)]) + #cbPalette[c(1:5,5,6,6)]) +
  #use free scales for each panel
  facet_wrap(~Scenario, nrow=2, scales='fixed', labeller= label_parsed) +
  scale_x_continuous(name=NULL, limits = c(2020,2099), breaks =c(2025,2050,2075,2099)) +
  scale_y_continuous("Change in global yield (%)", limits = c(-20,20))+#, breaks = c(-10,0,10,20)) + #, sec.axis = dup_axis())
  guides(colour = guide_legend(keywidth =1,keyheight =0.8, nrow=1, byrow=T, label.position = "right", label.hjust = 0, label.vjust = 0,
                               override.aes = list(size = override.size, linetype = override.linetype ))) #merge the legends for color and size
p0 <- p0 + theme_bw() + #theme_minimal() + #
  ggtitle("") + theme(plot.title = element_text(size=14, face="bold", hjust = 0.5)) +
  theme(panel.spacing.x = unit(1, "lines")) +  #facet spacing
  theme(strip.text = element_text(size = 13)) +
  theme(legend.text = element_text(size=12), legend.title = element_blank(), legend.position = "bottom") +
  theme(axis.text = element_text(size=13), axis.title=element_text(size=13))
# set transparency for both plot and panel background
p0 + theme(strip.background = element_blank(),
  panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
  panel.background = element_rect(fill = "transparent",colour = NA),
  plot.background = element_rect(fill = "transparent",colour = NA))

ggsave(paste0(wd1,"/figures/main/Fig3a-d.svg"), device ="svg",  #vertical plots
       units = "in", width = 7, height = 6.5, dpi = 600)



##save as PDF (1-column, 8.8cm)--resize line thickness and text sizes
casename = c("'SAI - RCP85'", "'MSB - RCP85'", "'CCT - RCP85'", "'Emissions reduction'")
casename = c(expression(''*bold(a)*'            SAI - RCP85'), expression(''*bold(b)*'            MSB - RCP85'), expression(''*bold(c)*'            CCT - RCP85'), 
            expression(''*bold(d)*'        Emissions reduction')) #expression('('*bold(i)*') RCP45 - RCP85'))
size <- c(1,1,1, 1, 1, 1, 1, 1)*0.4
override.size <- 1.1*0.4
p0 <-  eff.tot1.t[Variable%in%vars] %>%
  mutate(Var = factor(Variable, levels=vars)) %>%
  mutate(Scenario = factor(Scenario, levels=scenario, labels=casename)) %>%
  ggplot(aes(Year, Mean, group=Variable)) +
  #ggplot(aes(Year, Median, group=Var)) + #Weighted median shows better match with modeled yield change than weighted mean
  geom_line(aes(color=Var, linetype=Var, size=Var)) + #use Case as group and set override values in manual scale
  geom_hline(aes(yintercept=0), color="gray", size=0.3) +
  geom_linerange(aes(x=Year, ymin=Conf.L, ymax=Conf.U, color=Var), size=0.3, alpha=0.3, show.legend=FALSE) +
  #change linetype, color, size mannually
  scale_linetype_manual(values=linetype, guide = FALSE)+
  scale_size_manual(values=size, guide = FALSE) +
  scale_colour_manual(name="", labels=varname, values=varcolor)+ # cbPalette[c(1:4,4,5,5)]) + #cbPalette[c(1:5,5,6,6)]) +
  #use free scales for each panel
  facet_wrap(~Scenario, nrow=2, scales='fixed', labeller= label_parsed) +
  scale_x_continuous(name=NULL, limits = c(2020,2099), breaks =c(2025,2050,2075,2099)) +
  scale_y_continuous("Change in global yield (%)", limits = c(-20,20))+#, breaks = c(-10,0,10,20)) + #, sec.axis = dup_axis())
  guides(colour = guide_legend(keywidth =0.4,keyheight =0.3, nrow=1, byrow=T, label.position = "right", label.hjust = 0, label.vjust = 0,
                               override.aes = list(size = override.size, linetype = override.linetype ))) #merge the legends for color and size
p0 <- p0 + theme_bw(base_size = 7) + #theme_minimal() + #
  #ggtitle("") + theme(plot.title = element_text(size=10, face="bold", hjust = 0.5)) +
  theme(element_line(size = 0.2)) +
  theme(panel.spacing.x = unit(1, "lines")) +  #facet spacing
  theme(strip.text = element_text(size = 7, hjust=0.5)) +
  theme(legend.text = element_text(size=6), legend.title = element_blank(), legend.background = element_blank(), 
        legend.position = "bottom") + #legend.position = c(0.5, -0.15)) + #  
  theme(axis.text = element_text(size=7), axis.title=element_text(size=7)) + 
  theme(plot.margin=unit(c(0.1,0.2,-0.2,0.1),"cm")) #margins for top, then right, bottom and left
# set transparency for both plot and panel background
p0 + theme(strip.background = element_blank(),
           panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
           panel.background = element_rect(fill = "transparent",colour = NA),
           plot.background = element_rect(fill = "transparent",colour = NA))

ggsave(paste0(wd1,"/figures/main/Figure3a-d.svg"), device ="svg", units = "in", width = 3.5, height = 3)
ggsave(paste0(wd1,"/figures/main/Figure3a-d.pdf"), device ="pdf", units = "in", width = 3.5, height = 3)





#==================rescale/normalize all effects by radiative forcing (TOA radiative flux imbalance) to match that of ER
TOM <- read.table(paste0(wd1,"/data/NCL/TOM_radiative_imbalance.2020-2099.txt"), header = TRUE, sep=",")
TOM <- as.data.table(TOM)
TOM.ma <- TOM[, lapply(.SD, ma, order=5)] #5 year running mean as above
TOM.ma <- TOM.ma[, lapply(.SD, as.numeric)]
TOM.ma$Year <- 2020:2099
TOM.f <- TOM.ma[,.(ER=RCP45-RCP85, SAI=RCP85_SAI-RCP85, MSB=RCP85_MSB-RCP85, CCT=RCP85_CCT-RCP85), by=.(Year)]
#ER,SAI,MSB all have negative values, only CCT has a few positive values (RTMT > RCP85)
#original TOM radiative flux imbalance (not the same as sustained radiative forcing in equilibrum)
plot(ER~Year, data=TOM.f, type="l", col= "#26828EFF", ylim=c(-2,0.2))
lines(SAI~Year, TOM.f, col="#CC6600")
lines(MSB~Year, TOM.f, col="#56B4E9")
lines(CCT~Year, TOM.f, col="#B4DE2CFF")

#==Method: rescale the net effects by matching the mean and trend of SAI/MSB/CCT to those of ER
TOM.t <- TOM.f[, data.table(t(.SD), keep.rownames=TRUE), by=.(Year)]
names(TOM.t)[2:3] <- c("Scenario", "Data")

TOM.lm <- TOM.t[, .(intercept=lm(Data~Year)$coefficients[1], slope=lm(Data~Year)$coefficients[2],
                    Mean=mean(Data, na.rm = T)), by =.(Scenario)]
TOM.dtr <- TOM.t[TOM.lm, .(Data.dtr1=Data-i.slope*Year-i.intercept+i.Mean, #only remove the trend but keep the original mean
                           Data.dtr2=Data-i.slope*Year-i.intercept, #same as detrend() function, mean becomes zero
                           Data, intercept=i.intercept,slope=i.slope,Mean=i.Mean, Year), by=.EACHI, on =.(Scenario)]
lines(Data.dtr1~Year, TOM.dtr[Scenario=="ER"], col="#26828EFF", lty=2)
#lines(Data.dtr2~Year, TOM.dtr[Scenario=="ER"], col="#26828EFF", lty=2)
lines(Data.dtr1~Year, TOM.dtr[Scenario=="SAI"], col="#CC6600", lty=2)
lines(Data.dtr1~Year, TOM.dtr[Scenario=="MSB"], col="#56B4E9", lty=2)
lines(Data.dtr1~Year, TOM.dtr[Scenario=="CCT"], col="#B4DE2CFF", lty=2)

#adjust the mean and the rescale the slope of SAI/MSB/CCT to those of ER
ER.lm <- TOM.lm[Scenario=="ER",] #baseline
TOM.scaler <- TOM.lm[, .(mean.adj = ER.lm$Mean-Mean, slope.adj = ER.lm$slope-slope, 
                         mean.ratio = ER.lm$Mean/Mean, slope.ratio = ER.lm$slope/slope), by=.(Scenario)]
#scaling method 1: return detrended data to scaled values
TOM.scaled1 <- TOM.dtr[TOM.scaler, Data.new1:=Data.dtr2+slope*Year*slope.ratio + ER.lm$intercept, #add the intercept of ER
                      by=.EACHI, on=.(Scenario)]
#scaling method 2: avoid using intercept of lm, only need Year and original Data, mean(Year)=2059.5 (2020:2099)
TOM.scaled2 <- TOM.dtr[TOM.scaler, Data.new2:=Data+mean.adj+(Year-mean(Year))*slope.adj, #same as above method
                      by=.EACHI, on=.(Scenario)]
#but mean.adj is specific to the TOM data, cannot be applied to the new data eff.tot

#scaling method 3: only need Year and original Data
TOM.scaled3 <- TOM.dtr[TOM.scaler, Data.new3:=Data+mean(Data,na.rm = T)*(mean.ratio-1)+(Year-mean(Year))*slope.adj, #same as above method
                      by=.EACHI, on=.(Scenario)]
#method 3 is the best way to use (do not need to do lm regression and can directly rescale the mean and slope)

#now all scenarios have the same mean and slope 
TOM.scaled1[, .(Mean=mean(Data.new1, na.rm=T), slope=lm(Data.new1~Year)$coefficients[2]), by=.(Scenario)]
TOM.scaled2[, .(Mean=mean(Data.new2, na.rm=T), slope=lm(Data.new2~Year)$coefficients[2]), by=.(Scenario)]
TOM.scaled3[, .(Mean=mean(Data.new3, na.rm=T), slope=lm(Data.new3~Year)$coefficients[2]), by=.(Scenario)]
#(mean.ratio-1) can be applied to the mean of new data (eff.tot), 
#but slope.adj cannot be directly applied to partial and total effects because the slope of each effect variable is different from that of TOA imbalance


#======Finally use the ratio between the fit line of SAI/MSB/CCT and the fit line of ER as single scaler on any effect variable
eff.scaler <- TOM.scaled3[2022<=Year & Year<=2097, 
                         .(Year, scaler= lm(Data.new3~Year)$fitted.values/lm(Data~Year)$fitted.values,
                           fit=lm(Data~Year)$fitted.values, fit.new=lm(Data.new3~Year)$fitted.values), by=.(Scenario)] 
#the "scaler" scales the RTMT(TOA imbalance) mean and slope (i.e. the fit line) of SGs to that of ER
#all scenarios now have the identical fit line (same fit.new)



#==============supplementary figures S10-S11============
#plot TOM radiative imbalance of five scenarios
case = c("RCP45","RCP85","RCP85_SAI","RCP85_MSB","RCP85_CCT")

svg(paste0(wd1,"/figures/supplementary/Fig.S10a.svg"), width =6, height = 4) 
par(mar=c(2, 5, 2, 2) + 0.1)
plot(RCP85~Year, data=TOM.ma, type="l", col= casecolor[1], cex.main=1.0, lwd=2, xlab="", 
     main="(a) TOA radiative flux imbalance", ylim=c(0.5,3.5), ylab=expression("W m"^-2))
lines(RCP45~Year, TOM.ma, col=casecolor[2], lwd=2)
lines(RCP85_SAI~Year, TOM.ma, col=casecolor[3], lwd=2)
lines(RCP85_MSB~Year, TOM.ma, col=casecolor[4], lwd=2)
lines(RCP85_CCT~Year, TOM.ma, col=casecolor[5], lwd=2)
legend("topleft", c("RCP85", "RCP45", "SAI", "MSB", "CCT"), 
       lty=1, lwd=2, cex=1, bty="n", col=casecolor)
dev.off()


svg(paste0(wd1,"/figures/supplementary/Fig.S10b.svg"), width =6, height = 4) 
par(mar=c(2, 5, 2, 2) + 0.1)
plot(ER~Year, data=TOM.f, type="l", col= casecolor[2], cex.main=1.0, lwd=2, xlab="", 
     main="(b) TOA radiative imbalance difference from RCP85", ylim=c(-2,0.2), ylab=expression("W m"^-2))
lines(SAI~Year, TOM.f, col=casecolor[3], lwd=2)
lines(MSB~Year, TOM.f, col=casecolor[4], lwd=2)
lines(CCT~Year, TOM.f, col=casecolor[5], lwd=2)
lines(Data.new2~Year, TOM.scaled3[Scenario=="SAI"], col=casecolor[3], lty=2, lwd=2)  
lines(Data.new2~Year, TOM.scaled3[Scenario=="MSB"], col=casecolor[4], lty=3, lwd=2)   
lines(Data.new2~Year, TOM.scaled3[Scenario=="CCT"], col=casecolor[5], lty=4, lwd=2)  

lines(fit.new~Year, eff.scaler[Scenario=="ER"], col=casecolor[2], lty=1, lwd=2) # "#26828EFF"
lines(fit~Year, eff.scaler[Scenario=="SAI"], col=casecolor[3], lty=1, lwd=1.5)  
lines(fit~Year, eff.scaler[Scenario=="MSB"], col=casecolor[4], lty=1, lwd=1.5)   
lines(fit~Year, eff.scaler[Scenario=="CCT"], col=casecolor[5], lty=1, lwd=1.5)  
legend("bottomleft", c("RCP45 - RCP85", "SAI - RCP85", "MSB - RCP85", "CCT - RCP85"), 
       lty=1, lwd=2, cex=1, bty="n", #col=c("#003366", "#CC6600", "#56B4E9", "#B4DE2CFF"))
       col=casecolor[-1])
dev.off()





#===apply scaler to eff.tot1  
eff.scaler[Scenario=="ER",Scenario:="RCP45 - RCP85"]
eff.scaler[Scenario=="SAI",Scenario:="SAI - RCP85"]
eff.scaler[Scenario=="MSB",Scenario:="MSB - RCP85"]
eff.scaler[Scenario=="CCT",Scenario:="CCT - RCP85"]
eff.tot1.scaled <- copy(eff.tot1.t)
eff.tot1.scaled <- eff.tot1.scaled[eff.scaler, scaler:=i.scaler, by=.EACHI, on=.(Scenario)]


#==apply the scaler to all partial and total effects including CLM modeled difference, but excluding base yield and crop area
eff.tot1.scaled <- eff.tot1.scaled[!Variable%in%c("yield1", "yield","yield2","area2","yield0","area"),
                                   .(Year, Mean=Mean*scaler, Conf.L=Conf.L*scaler, Conf.U=Conf.U*scaler),
                                   by=.(Scenario,Variable)]

#double check the scaling factor is consistent at any year for a given Scenario
tail(eff.scaler[Scenario=="CCT - RCP85"])
TOM.t[Year==2095 & Scenario=="ER",Data]/TOM.t[Year==2095 & Scenario=="CCT",Data] #1.48
TOM.scaled3[Year==2095 & Scenario=="CCT",Data.new3]/TOM.scaled3[Year==2095 & Scenario=="CCT",Data] #1.48
eff.scaler[Year==2095 & Scenario=="CCT - RCP85",fit.new]/eff.scaler[Year==2095 & Scenario=="CCT - RCP85",fit] #1.40
eff.tot1.scaled[Year==2095 & Scenario=="CCT - RCP85" & Variable=="Tot.eff",Mean]/eff.tot1.t[Year==2095 & Scenario=="CCT - RCP85" & Variable=="Tot.eff", Mean] #1.40
eff.tot1.scaled[Year==2095 & Scenario=="CCT - RCP85" & Variable=="T.eff",Mean]/eff.tot1.t[Year==2095 & Scenario=="CCT - RCP85" & Variable=="T.eff", Mean] #1.40

eff.tot1.new <- rbind(eff.tot1.scaled, eff.tot1.t[Variable%in%c("yield1", "yield","yield2","area2","yield0","area")])

casename = c(expression('('*bold(a)*') SAI - RCP85'), expression('('*bold(b)*') MSB - RCP85'), expression('('*bold(c)*') CCT - RCP85'),
             expression('('*bold(d)*') Emissions reduction')) #expression('('*bold(i)*') RCP45 - RCP85'))

p0 <-  eff.tot1.new[Variable%in%vars] %>%
  mutate(Var = factor(Variable, levels=vars)) %>%
  mutate(Scenario = factor(Scenario, levels=scenario, labels=casename)) %>%
  ggplot(aes(Year, Mean, group=Variable)) +
  #ggplot(aes(Year, Median, group=Var)) + #Weighted median shows better match with modeled yield change than weighted mean
  geom_line(aes(color=Var, linetype=Var, size=Var)) + #use Case as group and set override values in manual scale
  geom_hline(aes(yintercept=0), color="gray") +
  geom_linerange(aes(x=Year, ymin=Conf.L, ymax=Conf.U, color=Var), alpha=0.3, show.legend=FALSE) +
  #change linetype, color, size mannually
  scale_linetype_manual(values=linetype, guide = FALSE)+
  scale_size_manual(values=size, guide = FALSE) +
  scale_colour_manual(name="", labels=varname, values=varcolor)+ # cbPalette[c(1:4,4,5,5)]) + #cbPalette[c(1:5,5,6,6)]) +
  #use free scales for each panel
  facet_wrap(~Scenario, nrow=2, scales='fixed', labeller= label_parsed) +
  scale_x_continuous(name=NULL, limits = c(2020,2099), breaks =c(2025,2050,2075,2099)) +
  scale_y_continuous("Change in global yield (%)")+# limits = c(NA,40))+#, breaks = c(-10,0,10,20)) + #, sec.axis = dup_axis())
  guides(colour = guide_legend(keywidth =1,keyheight =0.8, nrow=1, byrow=T, label.position = "right", label.hjust = 0, label.vjust = 0,
                               override.aes = list(size = override.size, linetype = override.linetype ))) #merge the legends for color and size
p0 <- p0 + theme_bw() + #theme_minimal() + #
  ggtitle("") + theme(plot.title = element_text(size=14, face="bold", hjust = 0.5)) +
  theme(panel.spacing.x = unit(1, "lines")) +  #facet spacing
  theme(strip.text = element_text(size = 13)) +
  theme(legend.text = element_text(size=12), legend.title = element_blank(), legend.position = "bottom") +
  theme(axis.text = element_text(size=13), axis.title=element_text(size=13))
# set transparency for both plot and panel background
p0 + theme(strip.background = element_blank(),
           panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
           panel.background = element_rect(fill = "transparent",colour = NA),
           plot.background = element_rect(fill = "transparent",colour = NA))

#Fig.3a-d scaled by TOA radiative imbalance
ggsave(paste0(wd1,"/figures/supplementary/Fig.S11a-d.svg"), device ="svg",  #vertical plots
       units = "in", width = 7, height = 6.5)







# #=======If considering LUC: all partial effects of T,RD,RI,P,RH,CO2,LUC
# #=======aggregate log effects (weighted by area only) before translating to % effect give best results on total and partial effects (see eff.mc.wm.crop)
# #=======Do not use pct wm effects or resampling on pct effects, they will bias partial effects of RD and RI although total R.eff look ok, the sum of T.eff+R.eff+M.eff is also larger than Tot.eff!
# eff.tot4 <- copy(eff.tot.ma)
# eff.tot4 <- eff.tot4[Scenario=="RCP45 - RCP85", Tot.eff:=Tot.eff3]
# eff.tot4 <- eff.tot4[Scenario=="RCP45 - RCP85", Y.mod.dif:=Y.mod.dif3] 
# #==NOTE==above weighted mean of grid-level change from dif.8crop only include shared crop grids between SSP2 and SSP5, but not crops distributed in different locations
# #compared tot.eff and Y.mod.dif
# eff.tot4[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff), by=.(Scenario)]
# eff.tot4[Year >=2075 & Year<=2099 & Level=="Mean", mean(luc.eff), by=.(Scenario)]
# eff.tot4[Year >=2075 & Year<=2099 & Level=="Mean", mean(Y.mod.dif), by=.(Scenario)] #7-10% on yield
# 
# #==include all crop grid cells of SSP2 and SSP5 and calc relative change in their global average yield
# eff.tot4 <- eff.tot4[Scenario=="RCP45 - RCP85", Y.mod.dif:=(yield2-yield0)/yield0*100]
# eff.tot4 <- eff.tot4[Scenario!="RCP45 - RCP85", Y.mod.dif:=(yield-yield0)/yield0*100]
# #compared tot.eff and Y.mod.dif
# eff.tot4[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff), by=.(Scenario)]
# eff.tot4[Year >=2075 & Year<=2099 & Level=="Mean", mean(Y.mod.dif), by=.(Scenario)] #6-11% on yield
# 
# 
# #% effect on global yield
# eff.tot4.t <- eff.tot4[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Year,Scenario)]
# names(eff.tot4.t)[3:6] <- c("Variable", "Mean", "Conf.L", "Conf.U")
# #names(eff.tot4.t)[3:4] <- c("Variable", "Mean")
# 
# 
# #======ploting=======
# linetype <- c(1,1,1, 1,1,1,1,1, 1)
# size <- c(1,1,1, 1, 1, 1, 1, 1,1)*1
# override.linetype <- 1
# override.size <- 1.1
# casename = c(expression('('*bold(f)*') SAI - RCP85'), expression('('*bold(g)*') MSB - RCP85'), expression('('*bold(h)*') CCT - RCP85'), expression('('*bold(i)*') RCP45 - RCP85'))
# scenario = c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85")
# #casename = c("'SAI - RCP85'", "'MSB - RCP85'", "'CCT - RCP85'", "'Emissions reduction'")
# #for verifying the match between regression estimate and modeled change
# vars = c("T.eff", "R.eff", "M.eff", "co2.eff", "luc.eff", "Tot.eff", "Y.mod.dif") 
# varname = c("Temperature effect", "Radiation effect", "Moisture effect", expression('CO'['2']*' effect'), "LUC effect", "Total effect", "CLM modeled")
# #finally plot all partial effects
# vars = c("T.eff", "RD.eff","RI.eff", "P.eff", "RH.eff", "co2.eff", "luc.eff", "Tot.eff")#, "Y.mod.dif")
# varname = c("T effect", "RD effect", "RI effect", "P effect", "RH effect", expression('CO'['2']*' effect'), "LUC effect", "Total effect")
# 
# 
# p0 <-  eff.tot4.t[Variable%in%vars] %>%
#   mutate(Var = factor(Variable, levels=vars)) %>%
#   mutate(Scenario = factor(Scenario, levels=scenario, labels=casename)) %>%
#   ggplot(aes(Year, Mean, group=Variable)) +
#   #ggplot(aes(Year, Median, group=Var)) + #Weighted median shows better match with modeled yield change than weighted mean
#   geom_line(aes(color=Var, linetype=Var, size=Var)) + #use Case as group and set override values in manual scale
#   geom_linerange(aes(x=Year, ymin=Conf.L, ymax=Conf.U, color=Var), alpha=0.3, show.legend=FALSE) +
#   #change linetype, color, size mannually
#   scale_linetype_manual(values=linetype, guide = FALSE)+
#   scale_size_manual(values=size, guide = FALSE) +
#   scale_colour_manual(name="", labels=varname, values=varcolor)+ # cbPalette[c(1:4,4,5,5)]) + #cbPalette[c(1:5,5,6,6)]) +
#   #use free scales for each panel
#   facet_wrap(~Scenario, ncol=1, scales='fixed', labeller= label_parsed) +
#   scale_x_continuous(name=NULL, limits = c(2020,2099), breaks =c(2025,2050,2075)) +
#   scale_y_continuous("Change in global yield (%)")+# limits = c(NA,40))+#, breaks = c(-10,0,10,20)) + #, sec.axis = dup_axis())
#   guides(colour = guide_legend(keywidth =1, ncol=2, byrow=F, label.position = "right", label.hjust = 0,
#                                override.aes = list(size = override.size, linetype = override.linetype ))) #merge the legends for color and size
# p0 <- p0 + theme_minimal() + #theme_bw()
#   ggtitle("") + theme(plot.title = element_text(size=14, face="bold", hjust = 0.5)) +
#   theme(strip.text = element_text(size = 13)) +
#   theme(legend.text = element_text(size=12), legend.title = element_blank(), legend.position = "bottom") +
#   theme(axis.text = element_text(size=13), axis.title=element_text(size=13))
# # set transparency for both plot and panel background
# p0 + theme(
#   #panel.grid.major = element_blank(),
#   panel.grid.minor = element_blank(),
#   panel.background = element_rect(fill = "transparent",colour = NA),
#   plot.background = element_rect(fill = "transparent",colour = NA)
# )
# ggsave(paste0(wd,"/crop/supplementary/Fig1f.Alleffects-Bootstrap-log-wtm-MLR-5yrma-final.jpg"),  #vertical plots
#        bg = "transparent", units = "in", width = 3.5, height = 11, dpi = 600)








